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ABSTRACT 

We present a general framework for matching the point-spread function (PSF), photometric 
scaling, and sky background between two images, a subject which is commonly referred to 
as difference image analysis (DIA). We introduce the new concept of a spatially varying pho- 
tometric scale factor which will be important for DIA applied to wide-field imaging data in 
order to adapt to transparency and airmass variations across the field-of-view. Furthermore, 
we demonstrate how to separately control the degree of spatial variation of each kernel ba- 
sis function, the photometric scale factor, and the differential sky background. We discuss 
the common choices for kernel basis functions within our framework, and we introduce the 
mixed-resolution delta basis functions to address the problem of the size of the least-squares 
problem to be solved when using delta basis functions. We validate and demonstrate our al- 
gorithm on simulated and real data. We also describe a number of useful optimisations that 
may be capitalised on during the construction of the least-squares matrix and which have not 
been reported previously. We pay special attention to presenting a clear notation for the DIA 
equations which are set out in a way that will hopefully encourage developers to tackle the 
implementation of DIA software. 
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1 INTRODUCTION 

Difference image analysis (DIA) aims to measure changes, from 
one image to another, in the objects that make up a scene. In as- 
tronomy, the objects are typically point sources changing in bright- 
ness or moving on the sky. Astronomical images are formed on a 
discrete detector array, after the sky scene suffers attenuation, geo- 
metrical distortion and blurring by the atmosphere and optics, su- 
perimposed on a sky background, and corrupted by detector noise. 
All of these effects are to different degrees non-uniform across the 
scene and variable on a variety of timescales. Furthermore, pairs of 
images of the same scene may suffer small misalignments in posi- 
tion or scale, or gross rotational misalignments. 
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The changes in object properties that we wish to measure 
are thus entangled with changes in the sky-to-detector, or scene- 
to-image, transformation. A residual difference image, formed by 
simple subtraction of one image from another, is generally domi- 
nated by changes in the transformation. To extract the astronomical 
information, we must accurately model the changes in astrometry, 
throughput, background, and blurring between the two images. We 
may then make corrections to match these effects from one im- 
age to another and subtract to form "cleaner" difference images, or 
we may model the original images including changes in both ob- 
ject properties and image transformations. While current DIA tech- 
niques are based on the former approach, we advocate the latter. 

The model adopted to represent changes in the scene-to-image 
transformation must include the following differential (or correc- 
tive) components: 
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• A coordinate transformation between the coordinate systems 
of each image to correct for image misalignments and/or differ- 
ences in distortion. 

• A photometric scaling that corrects for the changes in the at- 
tenuating effects of the atmosphere (and possibly the telescope op- 
tics) and differences in exposure time. 

• A background offset that corrects for changes in the sky back- 
ground emission. 

• A convolution transformation that corrects for the changes 
in the image point-spread function (PSF) as a result of changes 
in atmospheric conditions and/or the telescope optics (e.g. focus 
changes). 

Note that all of these components model differential corrections, 
not absolute values (e.g. the convolution transformation models the 
change in PSF shape between images, not the PSF itself). 

The state of the art in DIA includes the components for pho- 
tometric scaling, sky background offsets, and PSF convolution in 
the DIA modelling process. Recent developments ( lBramichll2008l . 
from now on BOS) also include fractional pixel translations in the 
model. Other image misalignments (rotation, scale, shear, distor- 
tion) are corrected by pre-registering the images before application 
of DIA, usually involving image resampling. 

The framework for the cur rent approach to DIA was intro- 
duced bv lAlard & Luptoiil i ll 9981) (from now on A98) for matching 
a reference image to a target image. The convolution kernel (includ- 
ing the photometric scaling) to be applied to the reference image is 
decomposed into a set of basis functions, and the differential back- 
ground offset is included as a polynomial of the image coordinates, 
which converts the problem of finding the corrective components 
to a standard li near least-squares formulation. A follow-up paper 
bv lAlard 120001) (from now on AOO) showed how the spatial varia- 
tion of the convolution kernel can be modelled by multiplying the 
kernel basis functions by polynomials of the image coordinates. 
The kernel basis functions chosen by A98 and AOO are Gaussians 
of different widths, modified by polynomials of the kernel coordi- 
nates. The user must specify the number of Gaussian basis func- 
tions to be employed, their associated widths, and the degrees of 
the modifying polynomials. However, the optimal choice of pa- 
rameters for generating the kernel basis functions is not obvious, 
although some investigation int o this matter has been performed 
jlsrael. Hessman & Schui]|2007l) . 

It is clearly desirable to find a set of kernel basis functions that 
are inherently simple, thereby being specified by a minimal param- 
eter set, and yet that can model the kernel with sufficient flexibility. 
A step towards this paradigm was made by B08 with the proposed 
representation of the kernel as a discrete pixel array where the ker- 
nel pixel values are solved for directly. This approach limits the 
requirements on the user to specifying the kernel size (and shape), 
and the kernel model is maximally flexible in modelling the most 
complicated convolution kernels (e.g. telescope jumps). B08 show 
that the new formulation is capable of modelling fractional pixel 
translations as part of the convolution kernel, thereby relaxing the 
requirement on image registration such that images need only be 
aligned to the nearest pixel before application of DIA. Spatial vari- 
ation of the kernel is handled by interpolation of kernel and differ- 
ential background sol utions on a grid. 

Soon after BOS, iMiller. Pennvpacker & w"hi5 feOOSh (from 
now on M08) specified a set of kernel basis functions built from 
delta-functions centred at different kernel coordinates. This choice 
of basis functions leads to a solution that happens to be equivalent 
to the BOS solution (see Section [T2t . but it is specified such that it 



fits into the A9S framework of equations. MOS also included a poly- 
nomial spatial variation of th e delta-function coefficients to model 
the kernel spatial variation. lOuinn. Clocchiatti & Hamuvl | |2010|) 
"rediscovered" the MOS work, but failed to impose any control 
on the photometric scaling while also fixing the value of the cen- 
tral kernel pixel, leading to a sub-optimal kernel model that cannot 
freely model fractional pixel translations. 

The choice of kernel basis functions in the A98 framework is 
fully down to the developer/user. While the delta-function basis (or 
delta basis for short) is very compelling, the number of free param- 
eters grows quickly with the adopted kernel size. Hence it makes 
sense to choose some coarser functions in the outer p art of the ker- 
nel wh ere there is little variation or signal/amplitude. [Albrow et al.l 
( I2OO9I) introduce the idea of binned kernel pixels in the outer part of 
the kernel, which greatly re duces the number of kernel parameters, 
an d l Yuan & Akerlo3 bOOSi) introduce a bicubic S-splines basis. 

One of the assumptions in the A9S DIA framework is that the 
photometric scaling between the reference image and the target im- 
age is characterised by a single number, which may be a reasonable 
assumption for images covering a small field-of-view (FOV), where 
spatial variations in atmospheric transparency and airmass are gen- 
erally negligible. However, DIA is now being applied in projects 
that generate images coveri ng multiple squa re d egrees each (e.g . 
Palomar Tran sient Factory - iRau et alj|2009l and lLawetai]|2009l . 
PanSTARRS - iKaiser et alj20IOl) , where non-uniform transparency 
is common (due to passing clouds) and extinction varies from one 
edge of the image to another due to airmass gradients across the 
field. Extension of the DIA framework to a spatially varying pho- 
tometric scale factor is therefore a necessary generalisation in the 
application of DIA to these projects. 

In Section |2| we take the step of generalising DIA to be able 
to cope with a spatially varying photometric scale factor, while si- 
multaneously modelling the spatial variation of the kernel shape 
and differential background. In presenting this generalised formu- 
lation, we also take the opportunity to present a clear set of DIA 
equations, with user-friendly notation, grouped in a logical way. 
The original DIA formulations in the literature (A98; AOO) are not 
so transparent in this respect, and the MOS formulation where delta 
basis functions are introduced omits the consideration of pixel un- 
certainties, has difficult notation, and misses a number of important 
simplifications with respect to this kernel basis (see Section [T2l (. 
Discussion of the most popular choices for the kernel basis func- 
tions and their implications with regard to the DIA formulation is 
made in Section [3] where we also introduce the mixed-resolution 
delta basis functions. In Section|4| we validate our algorithm using 
simulated data and we demonstrate it using some real data. Sec- 
tion |5] has been written to provide some implementation and opti- 
misation hints for the DIA developer, and the methodology that we 
propose will help to make the DIA algorithms more feasible with 
respect to the increasing data volume (image sizes and numbers) 
from the latest generation of time-series imaging projects. Finally, 
we state our conclusions in Section|6| 



2 THE GENERAL DIFFERENCE IMAGE ANALYSIS 
FORMULATION AND SOLUTION 

In this Section, we derive a general theoretical formulation of the 
difference image analysis problem from which all previously pub- 
lished formulations arise as special cases. This generalisation al- 
lows us to exercise control separately over the spatial variation of 
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each kernel basis function, the photometric scale factor, and the 
differential sky background, as we show in Sections [2.1l & l273l 



2.1 Defining The Target Image Model 

We start as in BOS by considering a pair of registered images sam- 
pled on the same pixel grid, one being the reference image with 
pixel values Rjj, and the other the target image with pixel values Ijj, 
where / and j are pixel indices referring to the column ; and row j of 
the image. We denote the spatial coordinate system in these images 
by {x,y), and the (x, v) coordinates of the (;,y)th pixel by {xi,yj). 
Exact image registration is not strictly necessary, since the best for- 
mulations for the kernel model include corrections for translational 
(but not rotational or otherwise) image misalignments, which has 
the advantage of avoiding problematic image interpolation in many 
cases. 

As first formulated by A98, we construct the model M(x,y) for 
the target image as the reference image convolved with a spatially 
varying kernel K{u,v,x,y) (where u and v are kernel coordinates) 
plus a spatially varying differential background B{x,y): 



M{x,y) = [R»K]{x,y)+B(x,y) 



(1) 



We wish to determine the best-fit convolution kernel and differen- 
tial background, and to do this we must first make further assump- 
tions about their functional form. We note that since the reference 
image is part of the target image model, it may be desirable to also 
determine the reference image pixel values Rjj. However, finding a 
solution to this issue is outside the scope of this paper. 

A98 made the important step of decomposing the kernel into 
a set of basis functions thereby linearising the expression in Equa- 
tion [T] Subsequently, AOO generalised the kernel decomposition 
to include the spatial variation of the basis function coefficients, 
which facilitated the modelling of the spatial variation of the ker- 
nel. We form the same kernel decomposition: 



K{u,v,x,y) = aq{x,y) Kq(u,v) 
9=1 



(2) 



where Kg{u,v) is the qlh kernel basis function, aq[x,y) is the gth 
spatially variable coefficient, and is the number of kernel basis 
functions. 

A polynomial is a sensible choice of model for the spatial vari- 
ation of the kernel basis function coefficients since it respects the 
linearity of the decomposition in Equation|2l and by specifying the 
polynomial degree, one may control the amount of spatial varia- 
tion that is to be modelled. The polynomial form for aq{x,y) was 
adopted by AOO with the same degree for each basis function coef- 
ficient. We generalise this further by modelling each coefficient as 
a polynomial with individual degree dq, providing a flexibility that 
we require later on: 



aq{x,y)='^ aq,„„^{x)"' £,{y)" 

m=0 n=0 



(3) 



where the aq,„„ are polynomial coefficients for the q\h kernel basis 
function. The coordinates (77 (x) , (y)) are normalised spatial coor- 
dinates defined by: 



^{x)- 



■{x-Xc)/N, 
: {y-yc)/Ny 



(4) 
(5) 



which follow from the Taylor expansion of the spatial coordinates 
{x,y) around the image centre {xc,yc) for an image of size Nx x A',. 



pixels. This coordinate conversion improves the orthogonality of 
the spatial polynomial termj^, and it prevents the significant poly- 
nomial coefficients from becoming progressively smaller for the 
higher order polynomial terms. 

As in A98, we also adopt a polynomial model of degree for 
the differential background: 



k=0 1=0 



(6) 



where the bki are the polynomial coefficients. 

We now have a model M{x,y) for the target image that is a 
linear combination of functions of x and y. This is easily shown by 
substituting Equations |2l [3] &|6] into Equation [T] and using the fact 
that convolution is distributive: 

M(x,y)=Y^[R®Kq][x,y)Y^ £ 77 W" ^ (j)" + £ ^ hm{xf ^{y)' 
q=\ m=0 n=0 k=Q 1=0 

(7) 

The target image is a discrete image of pixel values lij and 
therefore we wish to evaluate the model for the target image at the 
discrete pixel coordinates {xj,yj). Let us use Mjj to represent the 
discrete model image M{xi,yj) and {rii,^j) to represent the dis- 
crete coordinate array (t](x,), (yy)). Then, using the fact that the 
convolution of the reference image with the continuous kernel 
basis function Kq{u,v) is equivalent to a discrete convolution (see 
Appendix A), we have: 



d^, di,—m 



dn d-a-k 



Mi. 



with: 



l[^®'^],vi L aq^„r^'r^j+L L bkinn'j (8) 



9=1 



m— n— 



^=0 1=0 



(9) 



where r and s are pixel indices corresponding to the column r and 
row s of the discrete kernel basis function Kqrs defined by: 



Kq(u,v) du dv 



(10) 



We refer to [R(X) Kq]ij as a basis image since it is the linear 
combination of these basis images modified by spatial polynomi- 
als and combined with the differential background that constitutes 
the target image model. A basis image [R® Kq]ij is calculated from 
the discrete convolution of the reference image Rjj with the corre- 



sponding discrete kernel basis function Kq, 



. Equation 121 which 



implies that the reference image Rjj must extend beyond the pixel 
domain of the target image I/j. The discrete kernel basis function 
Kqrs may be defined directly, or calculated by analytical or numer- 
ical integration of Equation [TO] given a definition for Kq{u,v). Note 
that the terms for modelling the differential background in Equa- 
tion [8] can be thought of as multiplying a basis image that is set to 
unity at all pixels. 

All that is now required to fully define the model for the target 
image is to make a choice of suitable kernel basis functions, from 
which the corresponding basis images are derived. This is where 

' Although not considered here, further orthogonalisation of the spatial 
polynomial terms could be achieved by using, for example, Gram-Schmidt 
orthogonalisation. However, the orthogonalisation can only ever be approx- 
imate as the dot products that define orthogonality use inverse-variance 
pixel weights, and the variances depend on the model being fitted (see Sec- 
tionlZSt. 
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different authors have made different choices (e.g. the Gaussian 
basis functions, the delta basis functions, etc.), and we leave the 
treatment of these choices to Section [3] where we consider their 
implications in more detail. 

2.2 The Kernel Model 

Assuming that we have a solution for the polynomial coefficients 
aqmn of the kernel basis functions, we would like to know how to 
construct the discrete kernel model K^sij at any pixel {i,j) in the 
target image. This is achieved by defining: 

Kr.uj = j I J ^ K{u,v,Xi,yj) diidv (11) 
which, on substitution of Equations |2l[3]&[To] reduces to: 



Krsij = ^ Kgi-s ^ ^ aqmn Tji" ^'j 
q—l m— «— 



(12) 



2.3 Controlling The Spatial Variation Of The Photometric 
Scale Factor 

The kernel sum P,y = Y,rs^rsij-, which in general is a function of 
spatial pixel j), defines the photometric scale factor helween the 
reference image and the target image: 



Pij = ^ XI H H aqmn Tji ^ j 



(13) 



'9=1 



m— n— 



Our current formulation of the DIA problem in Section lTTl is 
such that Pij will vary across the image as a polynomial of degree 
equal to the maximum of the set of degrees rf^ax = max^y {dq} for 
the coefficients of the (sub-)set of kernel basis functions that have a 
non-zero sum. This can be seen by swapping the summation order 
in Equationll3land combining the kernel basis function coefficients 
into a single set of coefficients a^„: 



dm-ix dm„-m 



m=0 «=0 



where: 



52 '^g""> H 



(14) 



(15) 



This behaviour may be undesirable if we wish to employ a 
different degree of spatial variation in the photometric scale factor 
to the degree of spatial variation of the shape of the convolution 
kernel. AOO noted that those kernel basis functions with zero sums 
do not contribute to the spatial variation of the photometric scale 
factor, regardless of the spatial variation of their coefficients, and 
that one may always construct a new set of kernel basis functions 
that are a linear combination of the original set of basis functions. 

We assume that our kernel basis functions have been nor- 
malised to a sum of unity, or have a zero sum, and that our first 
kernel basis function KTi,-,, without loss of generality, has a sum of 
unity. We then form a new set of kernel basis functions as follows: 



f^qrs 



if g = 1 or Y.rs Kirs = 

ifq> \ and Kqrs 







(16) 



It follows that all of our new kernel basis functions if^,., have zero 
sums except for the first basis function K^^^ which has a sum of 
unity. 



Adopting our new set of kernel basis functions and dropping 
the prime from our notation, the photometric scale factor Pij re- 
duces to: 



di d\—m 

Bi=0 /1=0 



(17) 



which is a polynomial in the spatial coordinates (x, y) of degree d\ . 

Hence, by transforming the kernel basis functions as outlined 
above, one may specify a polynomial degree d\ of spatial variation 
for the photometric scale factor, associated only with the coeffi- 
cient of the first kernel basis function, and which we redefine as 
the degree rfp. Collectively, the spatial variation of the kernel ba- 
sis functions describes the kernel shape variations, and therefore 
the polynomial degree of spatial variation for the kernel shape is 
set by the value of max^jt/fy}, which is always greater than or 
equal to rfp. This is an important point to understand since if one 
wants to model the situation where the kernel shape is expected 
to spatially vary with a smaller degree than the photometric scale 
factor, then one should still fit a model with min^y {dq^ = dp. For 
example, to model the situation where the kernel shape is spatially 
invariant between two images but the spatial transparency pattern 
varies linearly (e.g. because of changes in airmass gradient), then 
one must adopt a linear spatial variation for all of the kernel basis 
functions. This enables the spatial variations of the zero-sum ker- 
nel basis functions to offset the spatial variations in kernel shape 
induced by the spatial variations of the unit-sum kernel basis func- 
tion. 

To summarise, we have shown how to decouple the spatial 
variation of the photometric scale factor from the kernel shape vari- 
ations (with the aforementioned caveat), which leads to three natu- 
ral types of spatial variation in the DIA formulation; namely, pho- 
tometric scale factor variations, differential background variations, 
and kernel shape variations, characterised by the degrees dp, d^, 
and ds = maxq {dq} >dp, respectively. 

2.4 Fitting The Target Image Model 

In order to fit the model in Equation[8]to the target image, we con- 
struct the chi-squared: 



I 

•j 



(18) 



where the (Jij represent the target image pixel uncertainties. Min- 
imising the chi-squared in Equation [Ts] falls into the class of gen- 
eral linear least-squares problems, since the model in Equation[8]is 
linear with respect to the unknown coefficients aq,,,,, and Z?^; to be 
determined. This class of problems has a standard solution proce- 
dure by construction of the normal equations. We refe r the reader 
to the treatment of this subject in Numerical Recipes jPress et al] 
|2007|) for more details. 

The normal equations are most compactly represented by the 
matrix equation: 



Ha = /3 



(19) 



where the square matrix H is the least-squares matrix, the vector a, 
is the vector of model parameters, and /3 is another vector. 

For each kernel basis function, there are 
!^q = [dq + \ ){dq + 2) /2 polynomial coefficients a^,,,,,, and for the 
differential background, there are A'b = (ds + 1)('^b +2)/2 poly- 
nomial coefficients fe^r./, leading to a total of A'par = (T^q^q) 
parameters to be determined. Hence the least-squares matrix H is 
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of size Np.^y by A'pa, elements, and the vectors a and j3 are of length Mij by using Iij, which enables the calculation of the initial ker- 
nel and differential background solution. In subsequent iterations, 
the current image model defined by Equation [S] should be used to 
set the Oij as per Equation [23] In Appendix B, we use an example 
to demonstrate the bias that can be introduced into the model pa- 
rameters if the iterative fitting procedure is not performed (see also 
SectionliH. 

It is also desirable to employ a A:-sigma-clip algorithm in order 
to prevent outlier target image pixel values from influencing the so- 
lution, including those from variable objects and cosmic ray events. 
This may easily be achieved by calculating the normalised residu- 



A'par elements. 

If we take z as a generalised index for all of the free param- 
eters, then we are simply assigning a one-to-one correspondence 
/ : z •<-> {q,m,n,k,l) that specifies which coefficient, aq^n or bj^i, 
corresponds to the current element a~ of the vector of parameters 
a. This mapping may order the parameters in an arbitrary way, but 
the ordering is only important for the efficient computation of H 
and p if one does not pre-calculate all of the necessary polynomial 
and basis images (see Sections fS . 1 1 & \5?2\ . 

Following from the definition of the model for the target im- 
age in Equation [S] the elements of the least-squares matrix H (i.e. 
the coefficients in the normal equations) and vector /5 may now be 
written out explicitly in terms of the basis images: 



for a- = aqnin and a^- 



m n 



for a. 



^qmn 



and 0!j/ = hjt'i' 



Lijrif+"''^^+"'\R®Kq,yj/olj 



(20) 



for cUr = bj^i and (Xv 
for 0!-; = bi-i and a^i = bj^iii 



for a- = a, 



qmn 



for a. = bj^i 



(21) 



Cholesky factorisation of the symmetric and positive-definite 
matrix H, followed by forward and back substitution is the most 
efficient and numerically stable method (Golub & Van Loan 1 9961) 
for obtaining the solution a = S to the normal equations. Explicit 
calculation of the matrix inverse H^' is only strictly necessary if 
one requires the covariance matrix cov{(X-,a^i) = {H^'}_,. We 
note that the calculation of the uncertainties in the elements of CI is 
one such case since the uncertainty a~ in each Sj is given by: 



als Eij = {Iij — Mij)/ajj and ignoring any pixels with |ey| > A: in 
subsequent iterations. The reliability of the ^-sigma-clip algorithm 
depends heavily on the accuracy of the adopted noise model, and 
since the initial (T,j values are calculated using an approximation 
to Mjj, we recommend that the sigma-clipping commences at the 
second iteration. 

Our final note in this Section is that the noise model in Equa- 
tion [23] could be improved, specifically by considering the noise 
introduced by the reference image, which is non-negligible when 
the S/N of the reference image is similar to that of the target image. 
AOO and BOS have previously considered such a noise model. Here, 
we explicit a useful noise model for a target image and a combined 
reference image that have been registered to the nearest pixel (i.e. 
avoiding image resampling): 



with: 



F2 .. 



+ 



Ml. 



GK 



2 

rcf,(!+r)(j+.s) 



m2 E 



A',- 



(24) 



(25) 



where the iJ^.^y represent the N^^, images that have been combined 
to create the reference image, and F,a,..,j and F^^fj^jj are the master 
flat-field images corresponding to the target image and constituent 
images of the reference image, respectively. 



{H^'lzz (22) 



2.5 The Noise Model And Iteration 

The calculation of the least-squares matrix H and vector /5 requires 
the adoption of a suitable noise model for the target image pixel 
uncertainties CT,j. BOS specify one such model as: 



' GFu 



■ + 



(23) 



where cTq is the CCD readout noise (ADU), G is the CCD gain 
(e~/ADU), and F,y is the master flat-field image. This model as- 
sumes that both the master flat-field image Fij and the reference 
image Rij are noiseless, which is a reasonable assumption for such 
typically high signal-to-noise (S/N) images. 

Most importantly, we note that in this noise model, the un- 
certainties Gij depend on the target image model Mjj and conse- 
quently, fitting Mij as described in Section |2!4l becomes an itera- 
tive procesfl In the first iteration, it is appropriate to approximate 



2.6 The Input Data 

Ideally, every pixel in the target image should be used in the calcu- 
lation of H and )3 , and therefore contribute to the kernel and differ- 
ential background solution. However, due to the nature of the con- 
volution process, the target image model is undefined in a border of 
width equal to half the kernel width around the image edges if the 
reference image is the same size as the target image, and therefore 
these target image pixels cannot be used in the calculation of H 
and /S. Also, "bad" pixels (e.g. bad columns/rows, hot pixels, satu- 
rated pixels, cosmic-ray events, etc.) should be excluded from the 
calculations, which means that any target image pixel [i, j) to be in- 
cluded in the calculation of H and j3 should be "good" in the target 
image, and that all reference image pixels to be used for calculating 
the target image model at (;, j) should be "good" in the reference 
image. This implies that a bad pixel in the reference image can 
discount a set of pixels equal to the kernel area in the target im- 
age, and therefore, as suggested in BOS, bad pixels in the reference 
image should be kept to a minimum, and kernels with excessively 
large footprints should be avoided when there are bad pixels in the 



reference image (e.g. see Section 2.3 of lBramich et al.ll201 ih 



^ Strictly speaking, the fact that the uncertainties (J,, depend on the target 
image model M,y also implies that minimising is no longer equivalent to 
maximising the likelihood. The maximum likelihood estimator is obtained 



instead by minimising x + '"(ff,^), which renders the fitting of the tar- 
get image model as a non-lineai' problem. 
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The areas of the target image which contain only sky back- 
ground and no astronomical objects will only contribute informa- 
tion on the differential sky background coefficients in the target 
image model. Hence, one may limit the set of target image pixels 
to be used in the calculation of H and /3 to a set of image sub- 
regions encompassing the higher S/N objects in the target image, 
which speeds the computations (fewer pixel values to be included 
in the required summations) while sacrificing some information. 
We note that contrary to the statements of some authors (e.g. M08), 
these sub-regions need not be centred on isolated stars. In fact, sub- 
regions of crowded high S/N objects (PSF-like or not) are precisely 
the image regions that contain the most information on the convo- 
lution kernel and differential background, because each pixel con- 
tains PSF and background information at a high S/N ratio. 

2.7 Difference Images 

We briefly mention that the definition of a difference image Dij is: 



Di 



-hi 



-Mi. 



(26) 



This image of residuals consists of noise (mainly Poisson noise 
from photon counting) and any differential flux from objects that 
have varied in brightness and/or position compared to the epoch 
of the reference image, since constant sources are fully subtracted 
during the DIA process. However, if an inappropriate kernel and/or 
differential background model is chosen, then unwanted system- 
atic errors will leave signatures in the difference image as large- 
amplitude high-spatial-frequency residuals at the positions of the 
brighter objects (for inappropriate kernel models), and as lower- 
amplitude low-spatial-frequency deviations in the difference image 
background from zero (for inappropriate differential background 
models). We note that if a reliable noise model exists, then the nor- 
malised difference image e,j defined by: 

. _iu-Mu 



(27) 



acts as a useful guide to the level of flux variation in any one pixel, 
since the pixel values in this image are in units of sigma-deviations. 

The purpose of producing a difference image is to enable accu- 
rate differential photometry to be performed in the absence of PSF 
crowding for all objects of interest (constant and variable). The ob- 
ject positions are presumed known from analysis of the reference 
image or from fitting of the differential flux on the difference im- 
age. 



3 COMMON BASIS FUNCTION CHOICES 

In this Section, we elucidate the common choices for the kernel ba- 
sis functions. We stress that since the choice of basis functions is 
fully independent of the DIA framework presented in the previous 
Section, the generation of a set of basis functions may be imple- 
mented as code that is completely separate from the DIA code. 



3.1 The Gaussian Basis Functions 

A98 introduced the Gaussian basis functions as a set of two- 
dimensional radially-symmetric Gaussian functions of different 
widths, each one modified by a polynomial of the kernel coordi- 
nates of a certain degree. The justifications for this choice are that 
an instrumental PSF is approximated by a Gaussian to first order, 
the convolution of a Gaussian by a Gaussian is also a Gaussian, and 



that a Gaussian decays rapidly beyond a given distance. The user is 
required to specify the number of Gaussian functions A'uau, and then 
for each Gaussian function (indexed by A), the user must specify 
the width Og^^ and the degree of the modifying polynomial D^^^^ . 
It follows that the definition of the gth kernel basis function corre- 
sponding to the Ath Gaussian with a modifying polynomial term of 
de gree ^/irau.H and degree ^/gau,i' ii^ the u and v coordinates, respec- 
tively, is given by: 



(28) 



gau,A 



The number of kernel basis func- 



tions Nk in this prescription is given by: 



Nk-- 



■t 

A = l 



(29) 



The Gaussian basis functions need to be numerically inte- 
grated via Equation [To] to form the corresponding discrete kernel 
basis functions, and then subsequently they should be transformed 
as detailed in Section l23] to allow control over the spatial variation 
of the photometric scale factor. Finally, we note that the adoption 
of a set of Gaussian kernel basis functions does not provide any 
simplification in the calculation of the basis images [i?® Kq]ij via 
Equation |9] 

Typical specifications for the Gaussian basis functions in 
the literature usually include three Gaussian functions, and the 
ISIS2 . ^ software developed by A98 and AOO adopts Gaussian 
widths of 0.7, 2.0, and 4.0 pix with modifying polynomials of de- 
grees 6, 4, and 3, r espectively, by de fault, resulting in 53 Gaus- 
sian basis functions. Ilsrael. Hessman & Schuh (2007) investigated 
how the optimal choice of Gaussian basis functions depends on the 
properties of the images for which DIA is to be performed (e.g. 
seeing, S/N, etc.), and although they manage to give some general 
recommendati ons, there seems to be n o unique answer. It has also 
been noted bv lYuan & Akerloi l l2008i) that the radial symmetry of 
the Gaussian functions may not be appropriate for elliptical PSFs, 
although it would be trivial to expand the Gaussian basis func- 
tion definition in Equation|28]to include elliptical two-dimensional 
Gaussians with an arbitrary centre and axis orientation. 



3.2 The Delta Basis Functions 

Let us introduce the definition of the Kronecker delta- function & , 



I if/ = j 
if ./ 



(30) 



Let us also assume that there exists a one-to-one correspondence 
g : q-ir^ {lJ.,v) which associates the ^th kernel basis function with 
the discrete kernel pixel coordinates (/i, v) such that, without loss 
of generality, q = {fl, v) = (0,0). Then we may directly define 
the ^th discrete kernel basis function Kqrs by: 



for g = 1 
for g > I 



(31) 



where we have already included the transformation as detailed in 
Section l23l to allow control over the spatial variation of the pho- 
tometric scale factor. It is clear that when q = 1, Ki^s obtains the 
value of I at {r,s) = (0,0) and elsewhere, and that when q> I, 



http://www2.iap.fr/users/alai'd/package.html 
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Kqrs obtains the value of 1 at (r, i) = (/x, v), — 1 at {r,s) = (0,0), and 
elsewhere. Thus adds flux to the PSF core, and the other Kqrs 
subtract flux from the core and add it back at displaced locations. 

We refer to this set of kernel basis functions as the delta basis 
functions. The set of delta basis functions may be chosen to cover 
any discrete kernel domain (e.g. circular - BOS, square - M08, etc.) 
by defining the number of kernel basis functions A'^ and the map- 
ping g appropriately. 

The basis images corresponding to the delta basis functions 
have a conveniently simple form that may be derived by substitut- 
ing Equation [3T] into Equation |9] and including a product of delta- 
functions to combine the two cases into one expression: 



Mixed-Resolution Delta Bosis Functions 



(32) 



ii+^l)U+v) 

Hence, the first basis image is the reference image itself, and the re- 
maining basis images are each formed by shifting the reference im- 
age by the appropriate integer-pixel shift, and then subtracting the 
non-shifted reference image. This has important speed and mem- 
ory implications when implementing the calculation of the least- 
squares matrix and vector (see Section|5}. 

BOS introduced the idea of solving directly for the kernel pixel 
values K,s of a spatially invariant kernel. We note that if we take 
for all q as an alternative definition for the discrete 



kernel basis functions in Equation|3T] then the corresponding basis 
images are given by [i?® Kq]ij = R{i+fi)(j+v]- This definition ig- 
nores any control that we may wish to exercise over the photomet- 
ric scale factor, but this is not an issue when considering a spatially 
invariant kernel (as in BOS). Substitution of this new result for the 
basis images into Equations |20]& [21] and assuming that the kernel 
and differential background are spatially invariant, leads directly 
to the least-squares matrix and vector derived by BOS from their 
direct solution approach. Hence, adoption of the delta basis func- 
tions in the A9S DIA framework is equivalent to solving directly 
for the kernel pixel values. A similar line of argument extends this 
conclusion to spatially variable kernels. 

The delta basis functions require minimal information from 
the user about the kernel shape and size for their specification. 
However, the dependence of the optimal kernel shape and size on 
the reference and target image properties has not yet been inves- 
tigated, although it is clear that the greater the difference in PSF 
width between the images, the larger the size of the convolution 
kernel that is required to match the PSFs. 



3.3 The Mixed-Resolution Delta Basis Functions 

The number of delta basis functions, and hence the number of co- 
efficients aq„,„, grows as the number of kernel pixels, which in turn 
grows as the square of the kernel radius. Since the least-squares 
matrix is a square matrix of size A'p^. by A'p^^ elements, the number 
of elements to be calculated in the least-squares matrix grows as the 
kernel radius to the fourth power. Hence, the time taken to calculate 
the solution for the coefficients a^„„ and increases considerably 

when solving for larger kernels. 

To address this performance issue. lAlbrow et al. I( l2009h intro- 
duced the idea of "binned" kernel pixels in the outer part of the 
kernel on the assumption that the kernel shows slower variations of 
smaller amplitude beyond a certain radius. Specifically, they intro- 
duced 3x3 binned kernel pixels beyond a kernel radius of 7 pix to 
replace the single kernel pixels, which greatly reduces the number 
of parameters to be solved for while maintaining a sufficiently large 
kernel. For example, for a circular kernel of radius 13 pix, which 
fits in a square array of 27 by 27 pixels, there are 577 single kernel 
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Figure 1. The distribution of single (red) and 3x3 binned (green) kernel 
pixels for a circular kernel of radius 1 3 pix that uses 3x3 binned kernel 
pixels beyond a radius of 7 pix. 



pixels. Adopting the 3x3 binned kernel pixels beyond a radius of 
7 pix results in 233 parameters, of which 177 are single kernel pix- 
els and 56 are 3x3 binned kernel pixels. The number of elements to 
be calculated in the least-squares matrix is consequently reduced to 
'-^16% without compromising the extent of the kernel model. Fig- 
ure[T]shows the distribution of single (red) and 3x3 binned (green) 
kernel pixels for this example. 

We generalise the idea of a binned kernel pixel to that of an 
extended delta basis function defined by: 



(1/^pix.?) + (5^0 5v0 - 1) 5r0 5i0 for {r,s) e Sq 
{SfiQ 5v0 - 1) 5,0 5so for (r, s) ^ Sq 



(33) 



where Sq is the set of kernel pixels spanned by the extended delta 
basis function (of any shape and spatial distribution), and A'pix.c/ is 
the number of elements in Sq. Note that we have assumed that the 
one-to-one correspondence g : ^ o (/i, v) is defined with = 1 <^ 
{fl,v) = (0, 0). Again, we have forced the sum of the extended delta 
basis function to be unity if it is the first kernel basis function (q = 
1), and to be zero if it is not (^ > 1), in order to be able to exercise 
control over the spatial variation of the photometric scale factor. 

The basis image corresponding to the extended delta basis 
function defined in Equation [33] is easily derived by substitution 
into Equation|9l 



[R(SKq]ij- 



{r.s)eS,i 



-{S^,oSvo-i)Rij (34) 



Therefore, this basis image is formed by averaging A'pix.fy versions 
of the reference image with each one shifted by the appropriate 
integer-pixel shift, and then subtracting the non-shifted reference 
image (except if this is the first basis image). Again, this has impor- 
tant speed and memory implications when implementing the calcu- 
lation of the least-squares matrix and vector (see Section [5}. We 
note that the basis images corresponding to the 3x3 binned ker- 
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nel pixels from lAlbrow et al] ( |2009|) are formed from integer-pixel 
shifted versions of a box-car smoothed reference image. 

We refer to a set of basis functions as mixed-resolution if they 
include any combination of delta basis functions and extended delta 
basis functions, and we emphasise that extended delta basis func- 
tions need not be square and they may be of any shape (e.g. circles, 
rectangles, rings, arcs, etc.). We also note that the delta basis func- 
tion is a special case of the extended delta basis function, and that 
overlapping extended delta basis functions are acceptable in a set 
of kernel basis functions as long as none of the extended delta basis 
functions may be constructed as a linear combination of any of the 
other kernel basis functions. If this condition is not met, then the 
solution for the coefficients in the target image model is degener- 
ate. Finally, we mention that mixed-resolution delta basis functions 
have the potential to be used in kernels with an adaptive resolution, 
which is a subject that has not yet been investigated in terms of its 
application to DIA. 



4 VALIDATING AND DEMONSTRATING THE 
ALGORITHM 

So far we have only examined the theory of our general DIA formu- 
lation. We now proceed to validate the algorithm using simulated 
images. We also demonstrate the ability of the algorithm to correct 
for a spatially varying differential transparency across the image 
area using real data. 

4.1 Simulated Image Data 

Our first task is to check that the algorithm can recover the exact 
model coefficients used to generate a set of simulated image data 
without any artificial noise added to the pixel values. By doing this 
we are validating our DIA formulation by confirming that there are 
no degeneracies in the target image model that we did not foresee. 

We generate a reference image of size lOOOx 1000 pix with a 
constant sky level of 1000 ADU and with 5000 stars. The stars are 
generated using a Gaussian PSF with a full-width half-maximum 
(FWHM) of 4 pix, pixel coordinates drawn from a uniform distri- 
bution over the detector area, and log-fluxes drawn from a uniform 
distribution between 2 and 5 (i.e. stars have fluxes between 10^ and 
10^ ADU). The image parameters that we have chosen are actually 
not important, and the tests in the absence of artificial noise give 
the same results so long as there are at least a few stars spread out 
over the image. 

We then generate a set of target images from the reference 
image using Equation [8] for various sets of kernel basis functions 
(Gaussian, delta, and mixed-resolution) and values for the corre- 
sponding coefficients, and for all combinations of dp, d^, and rfs 
(defined in Section [23t taken from the set {0,1,2,3}. We find that 
when we fit each target image with the model used to generate it, 
we can recover the exact input values (to within numerical preci- 
sion) of the coefficients aq,n„ and bj^i in Equation [8] for all cases. 
Hence we confirm that our algorithm works and that there are no 
hidden degeneracies. 

Next we generate a set of target images from the reference 
image by convolving the reference image with a spatially varying 
kernel of polynomial degree with the kernel normalised to a unit 
sum at each pixel. Then we multiply the convolved reference im- 
age with a polynomial surface of degree rfp representing the photo- 
metric scale factor and we add a polynomial surface of degree dg 
representing the differential background. We have done this for all 



combinations of rfp, rfg, and d^ taken from the set {0, 1,2}. In this 
set up, the degree of spatial variation of the kernel shape is actually 
df + since the polynomial surface for the photometric scale fac- 
tor multiplies the kernel pixel values which also spatially vary as a 
polynomial. Therefore, the appropriate (linear) target image model 
has rfp = dp, rfg =1^8' ''id rfs = d^+dg, and when we adopt such 
a model we find that we can recover the exact values for the model 
coefficients (again to within numerical precision). If we naively set 
dp = dp, d^ = dg, and d^ = d^ for our target image model, then the 
algorithm does not manage to perfectly fit the target image, leaving 
significant residuals. 

Now we consider how the algorithm performs for simulated 
images with added artificial noise. We adopt the same reference 
image as before and we use delta basis functions with di = dp = I 
and dq = ds = 2 for all q > 1. We define the kernel model to be 
a square array of 7x7 pixels. The target image model coefficients 
are arbitrarily chosen and specifically we set ai„„, = {1.1,0.3,0.1} 
for (m,;j) = {(0,0), (1,0), (0, 1)}. Also, we define 4=0 and set 
boo = 100. We then use all of these definitions in Equation [8] to 
generate a noiseless target image 5, ;. 

From the noiseless target image 5/j, we generate 10 noisy 
versions. Each noisy target image Ijj is formed by generating a 
1000x1000 pixel image S/y of values drawn from a normal dis- 
tribution with zero mean and unit o", and then computing: 

IiJ = Su+^ij^oi+Sij (35) 

where the coefficient of is derived from Equation |23] for G = 
1 e^/ADU and Fij = 1. We adopt a reasonable value for the read- 
out noise of cTq = 5 ADU. For each noisy target image, we fit the 
same model used to generate the noiseless target image, employing 
the iterative scheme described in Section [231 (but without sigma- 
clipping). 

In the plots along the diagonal of Figure (2] we show the dis- 
tributions of the coefficients ai,„„ for (m,n) = {(0,0), (1,0), (0, 1)} 
and boo as derived from the fits to the 10"' noisy target images. 
The red and black histograms represent the coefficient distributions 
after the first and third iterations, respectively, and the reason for 
iterating the solution is clear; namely, approximating M,j with lij 
in the noise model in Equation [23] in the first iteration introduces 
a significant bias into the fitted coefficients (in this example boo is 
underestimated by ~1 ADU or ~1%; see also Appendix B). We 
also report in the plots the measured mean and standard deviation 
of each coefficient distribution after the third iteration. The mea- 
sured means of the coefficient distributions are an excellent match 
to the input coefficient values (no differences to at least 5 signifi- 
cant figures), and the measured standard deviations are an excellent 
match to the formal uncertainties in the coefficients reported by the 
algorithm (calculated via Equation [22] and displayed as "Sigma"). 
For the coefficient distributions after the third iteration, we fit a 
Gaussian with mean and sigma equal to the corresponding input 
coefficient value and the formal uncertainty in the coefficient, re- 
spectively, and we plot the Gaussian fits as the blue curves. One 
can see that the coefficient distributions follow the Gaussian distri- 
butions very well. 

In the off-diagonal plots of Figure[2| we show scatter plots for 
all of the coefficient pairs that can be formed from aioo, fi io> 'JlOl > 
and boo using the results of the fits to the lO-' noisy target images. 
The red and black points represent the fitted coefficients after the 
first and third iterations, respectively. In each plot we also display 
the formal 1 C7-error ellipses (blue curves) as provided by the covari- 
ance matrix of the fit (see Section [2!4t . It is encouraging to see that 
there are virtually no correlations between the target image model 
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Figure 2. Plots along the diagonal: Histograms of the coefficients ai,„„ for (m,n) = {(0,0), (1,0), (0, 1)} and fooo as derived from the fits to the lO' noisy 
target images. The red and black histograms represent the coefficient distributions after the first and third iterations, respectively. The blue curves are fitted 
Gaussian distributions centred on the input coefficient values and with widths equal to the formal uncertainties in the coefficients. Off-diagonal plots: Scatter 
plots for all of the coefficient pairs that can be formed from aioo. f no. "101. and boo using the results of the fits to the 10^^ noisy tai'get images. The red and 
black points represent the fitted coefficients after the first and third iterations, respectively. The blue curves are formal IcT-eiTor ellipses. 



coefficients aioo> '^llOi and ajoi associated with the spatial variation 
of the photometric scale factor, or between the differential back- 
ground coefficient fcoo and ano or aioi- Also, as expected, there is 
a strong anti-correlation between the zeroth-order coefficients for 
the photometric scale factor and the differential background, oioo 
and boo. 

The anti-correlation between fljoo and boQ is a well-known 
feature of DIA that occurs when the reference image includes a 
non-zero background level. The kernel basis functions with non- 



zero sums in the target image model (Equation [Sj serve to blur 
and scale the reference image, including the background level, and 
hence the model terms for the differential background must com- 
pensate for this effect in the opposite sense. The consequence is 
that if the photometric scale is overestimated, then the differential 
background will be underestimated to compensate, and vice versa. 
To minimise the amplitude of this anti-correlation, we recommend 
subtracting the sky background from the reference image before ap- 
plying DIA (as suggested by BOS in their Section 2.2), a procedure 
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which could also include the subtraction of the spatially varying 
components of the background (i.e. to flatten the background). 

The results of our investigations in this Section lead us to con- 
clude that our DIA algorithm is working exactly as expected for 
simulated images with added artificial noise. 

4.2 Real Image Data 

We demonstrate our DIA algorithm using a pair of calibrated 
images from a commercial telescope (Celestron 8-inch Schmidt- 
Cassegrain / = 2032 mm) and CCD camera (Kodak KAF- 
1603ME) with a pixel scale of --^1.8 arcsec/pix, a FOV of 
0.26x0.26 degrees, and no filter. We designate one of the images 
as the reference image and the other as the target image. Since both 
images are undersampled, it is necessary to pre-blur them before 
applying DIA. We blur the reference and target images with Gaus- 
sian convolution kernels of FWHMs 3.5 and 4.0 pix, respectively. 

The target image was chosen specifically because it was taken 
through light clouds. We cropped the images appropriately in or- 
der to register them to the nearest pixel. We display the reference 
and target images in the top two panels of Figure [3] with the same 
linear scale and dynamic range of 700 ADU. The black regions are 
masked pixels that cover saturated stars. To the right of each im- 
age, we show three magnified image stamps corresponding to the 
red boxes marked in each image. These image stamps contain some 
of the brightest stars in the images which are most suitable for in- 
specting the quality of the difference images in the following tests. 

We proceed to fit the target image using the reference image 
and a set of delta basis functions representing a square kernel ar- 
ray of size 9x9 pixels. After some experimentation with different 
values for dp, rfg, and ds, we find that the differential background 
is only satisfactorily modelled for rfg > 3. The resulting difference 
images for each combination of rfp and taken from the set {0, 1} 
and with rfa = 3 are displayed in the middle four panels of Fig- 
ure [3] all with the same linear scale. The complicated residuals in 
the differential sky background are apparent in all cases. 

For {df,d^,d^) = (0,3,0), the dominant residuals at the star 
positions show an under-subtraction of the star fluxes towards the 
top-right of the difference image, and an over-subtraction of the 
star fluxes towards the bottom-left, which is clearly due to the pres- 
ence of spatial transparency variations that are not modelled by the 
spatially invariant photometric scale factor. This is also the case for 
{dp,dg,ds) = (0,3, 1), but since the kernel model is allowed to vary 
in shape across the image area, the zero-sum delta basis functions 
try to mitigate the spatial transparency variations by moving flux 
from the reference image background to the star PSF for those stars 
whose fluxes are under-subtracted, and by moving flux from the 
star PSF to the reference image background for those stars whose 
fluxes are over-subtracted, resulting in smaller residuals at the star 
positions but with the residuals spread out over a larger area. This 
is most visible in the image stamps on the right which still show 
under- and over-subtraction of the star fluxes, but spread out over 
more pixels. Setting {dp,dB,ds) = (1,3,0) successfully removes 
the under- and over-subtraction of the star fluxes from the differ- 
ence images, but instead leaves positive-negative residuals at each 
star position whose orientation is a function of position, which is a 
consequence of not modelling spatial variations in the kernel shape. 

Adopting {dp,dB,ds) = (1,3,1) produces difference images 
where only the brightest stars can be seen to be mildly under- or 
over-subtracted, which is a much better result than what current 
DIA algorithms are capable of producing (i.e. the {df,dB,ds) = 
(0,3,1) case). It is quite possible that further improvements in 



the difference image residuals may be obtained by adopting even 
higher polynomial degrees for dp, d^, and ds, but a full optimisa- 
tion of the production of the difference image in our example is 
outside of the scope of this paper 

In the bottom two panels of Figure [3] we reproduce the fitted 
photometric scale factor and differential background as a function 
of position over the image area which show that the atmospheric 
transparency diminishes and the sky background brightens for the 
parts of the target image that are more affected by clouds. This re- 
sult is to be expected since clouds attenuate the incoming light from 
outside the Earth's atmosphere, but they also increase the local sky 
brightness by scattering ambient light (e.g. light pollution, moon 
light, etc.) back to the ground. 

However, to be absolutely sure that this observed anti- 
correlation is not an artefact of our modelling procedure, we per- 
formed the following test. We cut out ten well-distributed im- 
age stamps around bright stars from the reference image, and we 
also cut out the corresponding stamps from the target image. For 
each pair of image stamps, we proceeded to fit the target image 
stamp using the reference image stamp and the same kernel con- 
figuration that we used to model the full target image, and we 
adopted a spatially invariant kernel and differential background (i.e. 
[dp.dg.ds) = (0,0,0)). We compared the photometric scale factor 
and the differential background derived from each fit, which rep- 
resent robust local estimates of these quantities, to the predicted 
values of these quantities at the stamp coordinates from our model 
for the full target image, and we found a very good agreement (to 
within ~2-4 per cent). This confirms that the results from our new 
DIA algorithm are fully consistent with the results that can be ob- 
tained using current DIA algorithms. 

This real data example has served as a proof-of-concept where 
we have demonstrated that we can use our DIA algorithm to suc- 
cessfully model a spatially varying photometric scale factor We 
have also shown how the results of solving for a spatially invari- 
ant kernel and differential background for small image sub-regions 
(stamps) in different parts of the image can be used to perform con- 
sistency checks on the solution for the target image model from our 
DIA algorithm. 



5 IMPLEMENTATION HINTS AND OPTIMISATION 
TRICKS 

Producing difference images is a very computationally intensive 
task, especially when modelling a spatially varying kernel as we 
have described in Section [2] However, many applications of DIA 
require quick (within seconds or minutes) processing of the target 
images (e.g. robotic searches for anomalie s in microlensing ev ents 
towards the Galactic bulge - RoboNet-II - iTsapras et alj|2009l su- 
perno vae searches - Palomar Transient Factory - Gal- Yam et al.l 
I20TI 

etc.). Hence the optimisation of the calculations required to 
produce the difference images is an important aspect of DIA. In the 
following subsections, we describe some useful optmisation tricks 
that may be used to obtain some substantial improvements in speed, 
and that have the potential to rival bru te force DIA imple mentations 
on graphical processing units (GPUs: lFluke et alll201 ih . 



5.1 Memory Considerations 

Firstly we consider what information needs to be stored in com- 
puter memory to enable the efficient calculation of the least-squares 
matrix H and vector /3, and the target image model Mij. We limit 
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Photometric Scole Foctor Differential Sky Bockground 




Figure 3. Top: A pair of calibrated images wliere the target image (right) was taken through light clouds. Middle: Difference images for various target 
image models. The red boxes are displayed as magnified image stamps to the right of each difference image. Note that in the bottom-left hand comer of the 
upper image stamp, the positive-negative residuals are caused by a moving object. Bottom: The spatial dependence of the fitted photometric scale factor and 
differential background for the target image model with {dp,d^,d^) = (1,3, 1). See text for more details. 
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ourselves to considering arrays that are of the size of the target im- 
age, since the other information that needs to be stored in computer 
memory (e.g. the discrete Icernel basis functions ic^„) takes up neg- 
ligible space in comparison. 

For efficiency, we need to pre-calculate and store in computer 
memory those images that will be used more than once in the calcu- 
lation of the difference image. In the general case, these images are 
the target image the reference image Rjj, the inverse- variance 
image 1 /c?-, the Nk basis images [R ® Kjlij, and the A'poiy polyno- 
mial images of the spatial coordinates T]" E," for m^n> 1 . If rf„,ax 
is the maximum degree of the polynomial spatial variation of the 
kernel basis function coefficients and the differential background, 
then the maximum degree of the polynomial images of the spa- 
tial coordinates in the least-squares matrix H is 2dniax (see Equa- 
tionl20t. which implies that: 

Wpoiy = [(2cU. + l)(2rf„,„ +2)/2] - 1 =rf„,ax(2rf„ax + 3) (36) 

For the Gaussian basis functions, with the typical choice of 
53 such functions (see Section UTb . it is perfectly feasible to store 
all of the corresponding basis images in computer memory (e.g. 
53 floating point 2000x2000 pixel images take up ~831 Mb of 
memory using IDL). Furthermore, the spatial variation of the kernel 
basis function coefficients is not usually modelled with a higher 
degree polynomial than a cubic polynomial, and a cubic polynomial 
variation requires A'poiy = 27. Again, it is possible to store all of the 
required 1 + 1 + 1+ 53 + 27 = 83 images in computer memory (e.g. 
83 floating point 2000x2000 pixel images take up ~I280 Mb of 
memory using IDL). 

For the delta basis functions, a typical circular kernel of ra- 
dius 10 pix generates 349 basis images, which is a more prob- 
lematic number of images to store in the computer memory (es- 
pecially for a 32-bit machine). However, the basis images for 
the delta basis functions may be generated without performing a 
computationally-costly convolution by simply subtracting the ref- 
erence image from a shifted version of itself (see Equation I32t. 
Hence, if one is prepared to recalculate each basis image as needed, 
and assuming that up to 27 polynomial images are required, then 
only I + 1 + I + 1+27 = 31 images need to be stored in computer 
memory. Similarly, using the same approach for mixed-resolution 
delta basis functions with A'^s resolutions only requires the storage 
of A'rcs versions of the reference image, each one produced by con- 
volving the original reference image with a box-car of the shape of 
the relevant extended delta basis function. 



5.2 Calculating The Least-Squares Matrix 

By far, most of the arithmetic operations required to fit the target 
image model and produce a difference image are performed in the 
construction of the least-squares matrix H and vector /5. In fact, as- 
suming that the inverse-variance, basis, and polynomial images are 
pre-calculated, and that D is the polynomial degree of spatial vari- 
ation of each kernel basis function and the differential background, 
then there are N^,„ = (W^ + 1)(D+ 1)(D + 2)/2 coefficients to be 
determined, and brute force computation of H requires the calcula- 
tion of yVp.„. entries, where the vast majority of these entries require 
3A'pix multiplications and A'pj^ — 1 additions (note that A'pj^ is the 
number of target image pixels that are being modelled). Further- 
more, /3 requires the computation of another A'par entries, where 
again the vast majority of these entries require 3A'pix multiplications 
and A'pix — 1 additions. Hence, the number of arithmetic operations 
A'op for the brute force computation of H and jS , normalised by A'pix , 



is given by: 

A'op~4iVpa.(A'p. + l) (37) 

We have already mentioned in Section |Z61 that limiting the 
target image pixels to be used in calculating H and /5 to a set of 
suitable image sub-regions minimises the number of required arith- 
metic operations for minimal loss of precision on the coefficients 
in the target image model. This clearly follows from the discussion 
in the previous paragraph. 

We have also noted in Section l2!4l that H is symmetric. This 
means that in reality only Npa,{Npra + I)/2 entries in H need to be 
calculated, and that the number of arithmetic operations reduces to: 

A'op~2yVp„-(A'p„ + 3) (38) 

Now we consider the order in which we may efficiently calcu- 
late the entries of H and /3 , and since H has the much larger number 
of entries, our choice is driven by the structure of H. Note that in 
the following, we treat the differential background as having a cor- 
responding basis image set to unity at all pixels (see Section l2Tt . 
Inspection of Equation [20] for H reveals that one has the choice of 
either: 

(i) For each of the (A'poiy + 1) polynomial images, cycle through 
the {Nk + 1)^ pairs of basis images to calculate the corresponding 
{Nk + 1)^ and Nk + 1 entries in H and /3, respectively. 

(ii) For each of the {Nf: + l)'^ pairs of basis images, cycle 
through the (A'poiy + 1) polynomial images to calculate the corre- 
sponding [(O + I ) (£) + 2) /2] ^ and (£) + I ) (£) + 2) /2 entries in H 
and p, respectively. 

We note that to calculate each entry in H, a pair of basis im- 
ages needs to be multiplied together before performing the re- 
quired summation, whereas the polynomial images are already pre- 
calculated from the coordinate images in computer memory, and 
therefore option (ii) is the most efficient because it minimises the 
number of the image multiplications that are required. Furthermore, 
in the case of the delta basis functions, the basis images are calcu- 
lated as needed, and therefore option (ii) also minimises the number 
of times that each basis image must be calculated. 

Having justified the choice of option (ii) for the order in which 
we should calculate the entries of H and /5 , we adopt a correspond- 
ing parameter ordering in the parameter vector a that leads to the 
structure for H that we illustrate in Figure glfor A't = 15 (arti- 
ficially small for clarity) and D = 2. The matrix H is made up 
of {N)c + 1)^ square sub-matrices (top panel of Figure |4ll, where 
each sub-matrix corresponds to the product of a single basis im- 
age pair [/?® Kq]ij [i?® Kq:]ij. Furthermore, each square sub-matrix 
has [(D+ l)(£' + 2)/2]^ entries, where each entry corresponds to a 
different polynomial image. However, within a single sub-matrix, 
there are only A'poiy + 1 = (-D + 1)(2£)+ 1) independent entries 
(Equation |36l bottom panel of Figure |4l(. In our specific example 
for D = 2, there are 15 independent entries out of 36 entries in each 
sub-matrix (i.e. less than half of the entries need to be calculated). 

The discovery of this property of H is exceptionally important 
because it greatl y decreases the number of required calculations. 
Neither M08 nor lOuinn. Clocchiatti & Hamuvl ( 1201 Ol) mention this 
optimisation, and AOO claim that the full modelling of the spatial 
variation of the kernel "quickly becomes intractable", and that "or- 
der 3 requires roughly 100 times more calculations than a constant 
kernel solution". We find that capitalising on the pattern in the sub- 
matrices of H for a spatial variation of the kernel of degree 3, one 
would only require A'poiy + 1 = 28 times more calculations than for 
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Figure 4. Top: The full least-squares matrix H divided up into (N^ + 1)^ 
square sub-matrices, where each sub-matrix has [(D+ l)(D + 2)/2]" en- 
tries. For this example we have adopted an artificially small value of 
= 15 for clarity, and D = 2. Bottom: A magnified view of a single square 
sub-matrix. Each sub-matrix in H has the same structure. Entries in the sub- 
matrix that employ the same polynomial image in their calculation have the 
same background colour (except for the single entries con'esponding to the 
polynomial images 1, 7],^, and ^J*). The polynomial term marked in each 
sub-matrix entry indicates the degree in the spatial coordinates (a,)') of the 
polynomial image correspondi ng to that entry. 
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Figure 5. Plot of the ratio of the number of arithmetic operations required 
to calculate H and /5 for our optimised algorithm compared to the same 
quantity for the brute force computation (black), and for the brute force 
computation that capitalises on the symmetry in H (red), when we adopt a 
set of delta basis functions representing a circular kernel. These ratios ai'e 
plotted as a function of the kernel radius (pix) and for D =0, 1, 2, and 3. 



a constant kernel solution, which is a very significant improvement 
in the potential performance of the algorithm. 

We are now in a position to develop an optimised algorithm 
for computing H and )5. We propose the following procedure: 

(i) For each row of square sub-matrices in H, carry out steps 
(ii)-(vii), and then finish. 

(ii) Calculate \R® Kq\ij / ofj and /,j [R® Kq]ij / <yfj for the cur- 
rent row, which requires 2A'pix multiplications. 

(iii) For each sub-matrix in the current row that lies on the di- 
agonal or in the upper half of H, carry out steps (iv)-(v), and then 
move on to step (vi). 

(iv) Calculate [R® Kq]ij[R® Kqi\ij / afj for the current sub- 
matrix, which requires A'pix multiplications. 

(v) For each pre-calculated polynomial image, calculate the ex- 
pression 77™+'" <^j'+" [R® Kq]ij [R® Kqi]ij I afj, which requires 
A'pix multiplications (except for m + m' + « + n' =0) and A'pix — 1 
additions, and fill out the relevant entries of the current sub-matrix. 

(vi) Fill out the entries of the sub-matrices in the current row 
that lie in the lower half of H by using the fact that H is symmetric, 
which takes a negligible number of operations. 

(vii) For each relevant pre-calculated polynomial image, calcu- 
late the expression '?;" ^']hj ^® K^];; / C7,^-, which requires A'pix 
multiplications (except for m + n = 0) and A'pix — 1 additions, and 
fill out the corresponding entries in jS. 

We now attempt to estimate the number of arithmetic op- 
erations that are required to calculate H and /3 using our opti- 
mised algorithm. Observe that step (ii) is repeated A'^^ + 1 times, 
steps (iv) and (v) are each repeated (A'^- + 1)(A'k: + 2)/2 times 
of which step (v) requires ~ (2A'pix)A'poiy + A'pix arithmetic op- 
erations, and step (vii) is repeated A'^: + 1 times and requires 
~ (2A'pix) [(£)+ l)(D + 2)/2] - A'pix arithmetic operations. Using 
^poiy = /)(2D + 3), then we derive the number of arithmetic op- 
erations in our optimised algorithm, normalised by A'pix, to be: 

A'op ~ (A'^ + 1 ) [a'^ (O + 1 ) (2D + 1 ) + 5Z)2 + 9Z) + 5] (39) 

In Figure[5] for D =0, 1, 2, and 3, we plot in black the ratio of 
the expression in Equation[39]to the expression in Equation[37]as a 
function of the kernel radius (pix) for a set of delta basis functions 
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representing a circular kernel. We see that for typical kernel radii 
of ~8-12 pix, we expect that our optimised algorithm will reach 
an efficiency in the number of arithmetic operations of ~0.251, 
0.167, 0.104, and 0.070 compared to the brute force computation 
forZ) =0, 1, 2, and 3, respectively. Also, in Figure|5] forZ) =0, 1, 2, 
and 3, we plot in red the ratio of the expression in Equation|39|to the 
expression in Equation [38] as a function of the kernel radius (pix) 
for the same set of delta basis functions. We further conclude that 
our optimised algorithm will reach an efficiency in the number of 
arithmetic operations of ~0.501, 0.334, 0.209, and 0.140 compared 
to the brute force computation that capitalises on the symmetry in 
H for D =0, 1, 2, and 3, respectively. 



6 CONCLUSIONS 

The general framework presented in this paper treats the problem 
of matching the PSF, photometric scaling, and sky background be- 
tween two images, where each of these components varies as a 
polynomial of the spatial coordinates. Where this paper improves 
over previous works on DIA are as follows: 

• We demonstrate how to model a spatially varying photometric 
scale factor within our framework, which is a new concept that will 
be important for DIA applied to wide-field imaging data that may 
suffer transparency and airmass variations across the field-of-view. 

• We show how to decouple the spatial variation of each ker- 
nel basis function, the photometric scale factor, and the differential 
background from each other, which allows more control over the 
level of spatial variation of each component in the target image 
model. 

• In Section|2]we develop what we hope is a clear notation and 
logical order for the DIA equations and methodology aimed at aid- 
ing others in creating DIA software implementations. 

• We prove the equivalence of adopting delta basis functions for 
the kernel model and solving directly for the kernel pixel values 
(BOS). 

• We introduce the mixed-resolution delta basis functions with 
the aim of reducing the size of the least-squares problem to be 
solved when using delta basis functions, and we elucidate their 
properties and implications for DIA. 

• We present some important optimisations in the calculation of 
the least-squares matrix which lead to a reduction in the number of 
arithmetic operations that need to be performed for typical kernel 
radii of ~8-12 pix to ~16.7%, ~10.4%, and ~7.0% compared to 
the brute force computation for linear, quadratic, and cubic spatial 
variations, respectively, of the target image model. 
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APPENDIX A 

Here we show that the convolution of the reference image Rjj with 
a continuous kernel basis function Kq[u,v) may be calculated as a 
discrete convolution. 

Firstly, consider the definition of continuous convolution ap- 
plied to the convolution of the reference image: 

[R(S)Kg]{x,y) = j j R(x + u,y + v)Kq(ii,v)dudv (40) 

where R{x,y) is a continuous representation of the reference image. 

Over the area of one pixel with coordinates (x,, Vj), the value 
of the reference image is a constant, i.e. R(x,y) = Rij for 
Xi — 1/2 < X < X/ + 1/2 and yj —\/2<y<yj->r\/2, and therefore: 

rs+l, rr+\ 

[R(SKg]{xi,yj)=Y,R(i+r){j+s) , " / ^Kg{u,v)dudv (41) 

where r and j are integer indices varying over the domain where 
the kernel basis function achieves non-zero values. 

Adopting the notation [R® for the image [R® K^](x,-,jy), 
then we may write: 

[R ® Kg] ij = Y.^(i+r)(j+s) Kjrs (42) 
rs 

/•J+i /•)-+! 

Kgrs= I , I , Kg{u,v)dudv (43) 

where r and s now represent the pixel indices corresponding to the 
column r and row s of the discrete kernel basis function Kq^s. 

Hence, the image [R® Kq](x,y) = [R® Kq]ij, which we refer 
to as a basis image, may be calculated via the discrete convolution 
defined in Eauationl42l 



APPENDIX B 

We wish to briefly investigate the consequences of approximating 
Mjj with Ijj in the noise model in Equation|23]as opposed to iterat- 
ing the solution and using the current image model from Equation[8] 
to update the noise model at each iteration. For this purpose we use 
the software developed in B08 for the case of a kernel and differ- 
ential background that are both spatially invariant. 

We create a 205x205 pixel noiseless reference image Rjj by 
setting a constant sky level of 1000 ADU and adding in 100 ob- 
jects, each of flux 10^ ADU and with a two-dimensional Gaus- 
sian profile of FWHM 4 pix, at random spatial coordinates drawn 
from a uniform distribution across the image area. We also cre- 
ate a 201x201 pixel noiseless target image Sij by convolving the 
Rjj with a discrete 5x5 pixel kernel calculated via numerical inte- 
gration of Equation [To] for a two-dimensional Gaussian of FWHM 
2 pix centred at the kernel centre and normalised to a sum of unity. 

We then perform the following experiment, adopting reason- 
able values for the readout noise and gain of Oq = 5 ADU and 
G = 1 e^/ADU, respectively: 

(i) We generate a 201x201 pixel image Z,y of values drawn 
from a normal distribution with zero mean and unit <J, and we con- 
struct a noisy target image /,y via: 

'iJ = Su + ^ij^<yi+Sij (44) 

where the coefficient of is derived from Equation [23] for G = 
I e^/ADUandF, , = 1. 
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(ii) We solve for a kernel and differential background that are 
both spatially invariant to match the reference image R^j to the tar- 
get image For the kernel model, we adopt 25 delta basis func- 
tions covering a 5 x 5 pixel array to match the actual domain of the 
discrete pixel kernel used to generate 5,^ from Rij. For the target 
image noise model (T/y we use Equation|23]with M,y approximated 
hylij. 

(iii) We record the photometric scale factor P\ and differential 
background B\ of the solution obtained in step (ii). 

(iv) We iterate the solution for the spatially invariant kernel and 
differential background three times (sufficient for convergence), 
each time using the current image model Mij calculated via Equa- 
tion [8] to set the target image noise model Oij via Equation |23] 

(v) Again we record the photometric scale factor Po and dif- 
ferential background B2 of the solution obtained during the final 
iteration in step (iv). 

We repeat the above experiment 10^ times and calculate the 
mean and standard deviation of each of the quantities Pi, B\, P2, 
and B2. We find that (Pi) - 1 = 5.38 x 10"* ± 1.68 x 10"^ and 
(Bi ) = - 1 .0085 ± 0.0020 ADU, where the uncertainty in the mean 
is estimated from the standard deviation divided by V 10^. The cor- 
rect solution in our experiment should have a photometric scale 
factor of unity and a differential background of zero. Clearly, solv- 
ing the DIA problem using the data to estimate the uncertainties on 
the pixel values in the target image introduces a bias of --^ 1 ADU 
in the differential background (and a very slight bias in the photo- 
metric scale factor). Hence one cannot assume that the background 
in the difference images produced using this method is zero, and 
aperture photometry on such difference images should include the 
computation and subtraction of a local background, and PSF pho- 
tometry should include the local background as a parameter in the 
fit. The bias in the differential background solution, which corre- 
sponds to an underestimated sky background in the target image 
model, is easily explained by the fact that the background pixels 
in the target image that randomly have smaller values than the true 
sky background are given more weight (or smaller uncertainties) 
in the fit than those background pixels that randomly have larger 
values than the true sky background. 

For the case where we iteratively solve the DIA problem us- 
ing the current image model to determine the uncertainties on the 
target image pixel values at each iteration, we find that {P2) — 1 = 
1.98 X lO^^i 1.68 X 10"^ and (63) = -0.0031 ± 0.0020 ADU. 
Therefore, at the precision of our experiment (which is well be- 
yond the photometric precision typically obtained for real data), 
we conclude that there is no bias in the derived photometric scale 
factor or differential background for this method, which validates 
the iterative method presented in Section [231 

Finally we mention that even though we only report one par- 
ticular experiment in this Appendix, we actually performed a range 
of experiments on artificial noisy target images generated with dif- 
ferent set-ups (e.g. different convolution kernels) and we found 
similar results in all cases. 
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